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Abstract 

We describe very fast electron dynamics for a graphene nanoribbon driven by a control electro- 
magnetic field in the terahertz frequency regime. The mobility as a function of bias field has been 
found to possess a large threshold value when entering a nonlinear transport regime. This value de- 
pends on the lattice temperature, electron density, impurity scattering strength, nanoribbon width 
and correlation length for the line-edge roughness. An enhanced electron mobility beyond this 
threshold has been observed, which is related to the initially-heated electrons in high energy states 
with a larger group velocity. However, this mobility enhancement quickly reaches a maximum 
governed by the Fermi velocity in graphene and the dramatically increased phonon scattering. 
Super-linear and sub-linear temperature dependences of the mobility are seen in the linear and 
nonlinear transport regimes, which is attributed separately to the results of sweeping electrons 
from the right Fermi edge to the left one through elastic scattering and moving electrons from low- 
energy states to high-energy ones through field-induced electron heating. The threshold field is 
pushed up by a decreased correlation length in the high field regime, and is further accompanied by 
a reduced magnitude in the mobility enhancement. This implies an anomalous high-field increase 
of the line-edge roughness scattering with decreasing correlation length due to the occupation of 
high-energy states by field-induced electron heating. Additionally, a self-consistent device mod- 
eling has been proposed for graphene transistors under an optical modulation on its gate, which 
employs Boltzmann moment equations up to the third order for describing fast carrier dynamics 
and full wave electromagnetics coupled to the Boltzmann equation for describing spatial-temporal 
dependence of the total field. Finally, a detailed comparison of the derived Maxwell-Boltzmann 
moment equations in this paper with the well known Vlasov-Maxwell equations is also included. 
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I. INTRODUCTION 



The engineering achievement of isolating graphene sheets- - - from graphite has inspired 
many studies aimed at understanding basic underlying physics- 1 ^ as well as finding possible 
applications to carbon-based electronics-. The low-field linear transport of charge carriers 
in a graphene layer, has received a considerable amount of attention. - — Recent reports 
on the successful fabrication of ultra-fast graphene transistors— and photodetector— has 
further advanced this research frontier into the fields of electronics and optoelectronics. 
The graphene transistor was reported to be used as both an electrical modulator— with a 
frequency as high as ~ 10 GHz and as a sensitive photo-detector for imaging—. However, 
similar investigations of linear transport in graphene nanoribbons (GNRs) have only been 
given relatively little attention so far.— ~— 

Early theoretical studies— 1 ^ on electron transport in graphene nanoribbons were re- 
stricted to the low-field regime, where a linearized Boltzmann equation was solved within a 
relaxation-time approximation. In this paper, the non-equilibrium distribution of electrons is 
calculated exactly by solving the Boltzmann transport equation beyond the relaxation-time 
approximation for nonlinear electron transport in semiconducting graphene nanoribbons. 
Enhanced electron mobility from initially-heated electrons in high energy states is antici- 
pated. An anomalous increase in the line-edge roughness scattering for large electric fields is 
obtained with decreasing roughness correlation length due to the occupation of high-energy 
states by field-induced electron heating. The semi-classical Boltzmann transport equation is 
expected to be applicable to the diffusive band-transport regime with relative smooth edges 
for graphene nanoribbons, instead of the hopping and tunneling between localized states 
with rough edges. 

Macroscopic simulation of semiconductor device physics has proceeded with solving cou- 
pled Maxwell-Boltzmann equations.— 1 ^ By employing the quasi-equilibrium Fermi-Dirac 
distribution in the Boltzmann transport equation, we can obtain the spatial dependence of 
both chemical potential and temperature, which paves the ground for drift- diffusion and 
hydrodynamic charge transport theories.— 1 ^ An ensemble of interacting electrons can be 
thermalized quickly with a high density, which justifies the assumption of a quasi-equilibrium 
Fermi-Dirac distribution for hot carriers. To determine the spatial dependence of chemical 
potential and temperature, device simulators usually couple Maxwell equations for the elec- 
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tromagnetic fields to conservation relations for carrier density, the current and energy. Sev- 
eral investigations on hydrodynamics are related to determining the proper mathematical 

representation of carrier thermal conductivity and the moments of the Boltzmann equa- 
te 22^25 

The outline of the remainder of this paper is as follows. In Sec. Ill A[ we solve exactly the 
semi-classical Boltzmann transport equation for low-temperature electron transport in semi- 
conducting graphene nanoribbons by including impurity, line-edge roughness and phonon 
scattering effects at a microscopic level. Based on the calculated non-equilibrium distribu- 
tion as a function of wave number along the ribbon, we present detailed numerical results 
for the electron mobility as a function of either the applied electric field or the lattice 
temperature for various impurity scattering strengths and correlation lengths for line-edge 
roughness. Our numerical results are presented in Sec. Ill Bl with some discussions. In Sec. 
IIII Al we derive the moment equations from the Boltzmann equation up to the third order 
for electron dynamics in n— doped graphene as a generalization to hydrodynamic model. At 
the same time, the self-consistent field equations are also derived in Sec. IIII Bl within the 
Maxwell-Boltzmann frame. Finally, the conclusions of this paper are presented in Sec. IIVI 

II. NONLINEAR TRANSPORT IN GRAPHENE NANORIBBONS 

In this section, we employ the Boltzmann transport model with inclusion of scattering at a 
microscopic level to study high-field nonlinear transport of electrons in graphene nanoribbons 
along with some numerical results. 

A. Nonlinear Boltzmann Transport Model 

Here, we investigate single subband nonlinear transport only in the armchair nanoribbon 
(ANR) configuration—. Low electron densities, moderate temperatures, ionized impurities 
and line-edge roughness are considered— >2L2&. As a result, negligible pair scattering—, 
optical and out-of-plane flexural phonons— , inter-valley scattering, volume-distributed and 
short-range impurity scattering— will all be neglected. Therefore, the electron-like branch 
for n— doped graphene can be represented on a k— space mesh,- via its dispersion and 
corresponding wave-function, as 
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e 5 = hv F x jk) + (ir/3Wy , (1) 

e i(2w/3a ~K)x ^ 



.1 

i 



Here, vp = 10 6 m/s is the Fermi velocity in graphene and L is the quantization length of the 
ribbon. For semiconducting ANR, k = ir/3W <C 27r/3ao is the quanta of the transverse wave 
vector and <f>j = tan -1 (kj/n) is the phase separation between the two graphene sublattices. 
The electron wave numbers kj = [j — (N + l)/2] 5k are given on the discrete mesh by 
j = 1, 2, . . . , N for large odd integer N, 5k = 2 fc max / (N — 1) is a small mesh spacing, and 
&max is chosen to ensure that scattering induced population of higher electron-like branches 
can be neglected. The minimum in the energy dispersion curve corresponds to the central 
mesh point j = M = (N + l)/2. Also, W = (Af + 1) a /2 is the width of a ribbon 
expressed in units of the size of graphene unit cell ao = 2.6 A and the number of carbon 
atoms M across the ribbon. According to the dispersion relation in Eq. (fTTh the electron 
group velocity Vj for semiconducting ANRs is given by Vj = Vp (Hi/pkj/Ej). Additionally, we 
assume that the electron-like branch is filled up to \kj\ = kp at zero temperature with the 
Fermi wave number and energy given by kp = 7mi£)/2 and £p = e(kp), respectively. For a 
chosen temperature T and chemical potential /io, the linear electron density in ANR follows 

JV 

from Hid = 5k/ir ^ /• , with = {1 + exp [(ej — fj, ) /ksT]}" 1 being the equilibrium 
Fermi-Dirac distribution function. 

Conventionally, the non-equilibrium carrier distribution function is partitioned as fj = 
fj ^ +gj. The deviation from the equilibrium Fermi distribution under a strong electric field 



is described by the set of reduced nonlinear Boltzmann transport equations 
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= evjfo l-^-j ~ 22 S jtj ,(t) gj ,{t) . (3) 

\ J J j'^M 

In this notation, g'At) = gj(t) — gj\i(t) is the reduced form of the dynamical non-equilibrium 

part— of the electron distribution function. The reduced form accounts for particle number 

conservation condition, i.e., 9j{t) — 0. Reduced scattering matrix elements «Sj = 

i=i 

~ Sj,M(t) are defined via its components 
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Sjj'it) = 5 hi' [ s j°j> + S j n fWf} + s f (! - ^-,(iv+i)/2)] - Sj +j > iN+ i [Sf (1 - 5 jt{N+ i )/2 )] 

Here, the elastic scattering rate is given by: 

5f = 5 - + = ( 70 + _2L_) ( - ) p + oos(2^)] , (5) 

where 70 ~ n2D denotes the impurity scattering rate at the Fermi edges and vf = Vj(kp)- 
Since the momentum difference between two valleys is rather large, only short-range impu- 
rities (such as topological defects) with a range smaller than the lattice constant will give 
rise to inter-valley scattering, ji = 2 {ixv F 5b/'5W 2 ) 2 (A /u F ) is the scattering rate due to 
edge roughens, with 5b ~ 5 A being the amplitude and A being the correlation length of 
the roughness. 

The dominating inelastic scattering mechanism is provided by longitudinal acoustic 
phonons at low temperatures. The static inelastic scattering rates are given by the fol- 
lowing matrix elements 



s fy = y k H s tr + > (6) 



2n 



Sf >f = 9{±eyTe j ) 



^al\ e o' ~ £ j 



2h 2 cjp M LWe^ F {\k j/ - k 



.1 1 



[l + cos(<^-<^)] . (7) 



Here, f~ = ff\ /+ = l-ff\ n j>f = N (\e y -e^/h), N (cu g ) = [exp(hw g /k B T) - 1]" 1 is the 
Bose-Einstein function for thermal equilibrium phonons; Dal ~ 16 eV is the deformation 
potential, pu ~ 7.6 x 10~ 8 g/cm 2 and c s ~2x 10 6 cm/s are the mass density and sound 
velocity in graphene. The scattering potentials are screened by free carriers and described by 
a dielectric function. Here, we assume that the inelastic scattering is shielded by the static 
Thomas- Fermi dielectric function in its general form— ^ eTF(\kf — kj\). The screening of 
elastic scattering potentials is given approximately by €tf ~ 1 + [e 2 l^ 2 eQe r hvp) under the 
metallic limit (2k F W > 1) with e r « 3.9. 

The nonlinear dynamical phonon scattering rate is 



> (8) 
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which is also responsible for the nonlinear electron transport and electron heating due to its 
dependence on g'y{t). 

Once the non-equilibrium part, g'j(t), of the total electron distribution function has been 
determined using Eq. fl3]), the transient drift velocity, v c (t), of the system can be calculated 
with the use of 



v c (t) 



N 



(o) 



E fe-^W) • ( 9 ) 

We note that the thermal-equilibrium part of the electron distribution does not contribute 
to the drift velocity. The steady-state drift velocity v& of electrons is given by v c (t) by taking 
the limit t — > oo. The corresponding steady-state conduction current is given by / = eninv^- 
The differential electron mobility for nonlinear transport is generalized to /i e = dv^/dJ^Q. 
Numerical simulation of these quantities in specific ANRs is presented below. 



B. Numerical Results and Discussion 



Figure DJa) presents our calculated electron mobilities function of applied electric 

field Jo at T = 10 K (blue curve) and T = 6K (red curve), respectively. Clearly, from 
Fig.[T](a), a strong /^-dependence for /i c appears at a lower value of jFq at higher tempera- 
ture T in the nanoribbon. This J-"o- dependent electron mobility fi c has its physical origins 
in the dynamical electron-phonon scattering rate <Sjy{g6} through its dependence on the 
electron distribution function given in Eq. (JSJ). Consequently, it is reasonable to expect a 
lower threshold field, J 7 *, for entering into a nonlinear transport regime (J 7 > J 7 *) due to 
enhanced nonlinear phonon scattering at T = 10 K. The value of J 7 * strongly depends on 
the parameters of the system, such as T, riin, 7o and A , and an analytic expression for J 7 * 
cannot be obtained in the nonlinear transport regime. On the other hand, as F$ — > 0, fi e 
is larger at T = 10 K than T = 6 K due to thermal population of high-energy states with a 
large electron group velocity. The initial decrease of /x e with T§ is attributed to the gradual 
increase of the frictional force from phonon scattering by jFq. At T = 10 K, /i e is roughly 
independent of Fq below 0.75kV/cm (linear regime). However, \x e increases significantly 
with jFq above 0.75kV/cm (nonlinear regime). Eventually, /i e decreases with J-~ once it 
exceeds 1.5kV/cm (heating regime), leading to a saturation of the electron drift velocity. 
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The electron group velocity \vj\ increases with the wave number for small \kj\ values, as can 
be seen from Eq. (pP). However, the increase of \vj\ slows down toward its upper limit vp 
provided \k^\ ^> ir/3W but still within the single-subband regime. The increase of /i e with 
J-"o i n the nonlinear regime comes from the initially-heated electrons in high energy states 
with a larger group velocity, while the successive decrease of /x e in the heating regime comes 
from the combination of the dramatically increased phonon scattering and the upper limit 
vp mentioned above. In Fig.[Tfb), the calculated electron drift velocities v& are plotted as a 
functions of temperature when T§ = 2kV/cm (blue curve) and F$ = lkV/cm (red curve). 
The fact that /x e increases with T monotonically in both cases implies the electron scattering 
in two samples is not dominated by phonons but by impurities and line-edge roughness. Dif- 
ferent behaviors in the increase of fi c with temperature can be seen from Fig.[TYb) for linear 
and nonlinear electron transport. At J-"q = 2kV/cm for the high-field nonlinear transport, 
t>d (or /i e ) increases with T sub-linearly. On the other hand, rises super-linearly with T 
for the low-field linear transport at = lkV/cm. These different T dependence in /x e for 
linear and nonlinear transports can be directly related to the non-equilibrium part of the 
electron distribution function gj. 

The effects due to impurity scattering are compared in Figs.[2]^a) and^b). In Fig.^a), 
we present a comparison of mobilities as a function of J-*o at T = 10 K for 70 = 1.0 x 10 13 s _1 
(blue curve) and 70 = 1.0 x 10 14 s~ 1 (red curve). As To — » 0, fi e is greatly reduced by 
strong impurity scattering with 70 = 1.0 x 10 14 s~ 1 in the linear regime. In addition, for 
70 = 1.0 x 10 14 s _1 , J 7 * is pushed upward from about 0.75kV/cm to 1.75kV/cm, leaving us 
with a roughly .^-independent fi e in this case for the whole field range shown in this figure. 
The comparison for T-dependence of [i e is presented in Fig.|2](b) for Tq = 2kV/cm, where 
a sub-linear increase of /i e with T for weak impurity scattering is switched to a super-linear 
relation in the strong impurity scattering case. 

Finally, the effect of correlation length for the line-edge roughness on the transport is 
demonstrated in Figs.[Ha) and[3](b) by fixing W = 50 A and changing Ao from 200 A to 50 A. 
As shown in Eq. (JH]), the line-edge roughness scattering can be either reduced or enhanced 
by decreasing Ao, depending on \kj\ 1/2 Ao or \kj\ 3> l/2Ao- For our chosen sample 
with riiD = 1.0 x 10 5 cm _1 , we find that the condition \kj\ 1/2A is satisfied in the 
low- field limit (\kj\ ~ kp), while \kj\ ^> 1/2Aq holds for the high field limit due to electron 
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heating. Therefore, we find from Fig.[3]^a) that fi c is increased as J-q — > when Ao is reduced 
to 50 A in the low field regime. However, the value of J 7 * for /i e is pushed upward for 
A = 50 A in the high-field regime, which is further accompanied by a reduced magnitude 
in the enhancement of /x e with J-" - This anomalous feature associated with reducing A 
also has a profound impact on the T-dependence of fi c as shown in Fig.[^b), where the 
increasing rate of /i e with T in the high field regime becomes much lower with Ao = 50 A 
than for A = 200 A. 



III. OPTICAL MODULATION TO GRAPHENE LAYERS 

In this section, we will first derive the Boltzmann moment equations up to third order 
in describing very fast electron dynamics driven by a control electromagnetic field in the 
terahertz frequency regime. At the same time, a self-consistent equation is also derived for 
the total driven field, including the induced optical coherence in a graphene layer. 



A. Boltzmann Moment Equations 

For electrons in a two-dimensional n— doped conducting graphene layer with its band 
structure given by a Dirac cone, the Boltzmann equation for the electron distribution func- 
tion k||, t) around K (or K') valley point is 



df(r h kn, t) dry 



■ Vk||/(rii, kn, t) = - 



(10) 



coll 



at dt "' ' dt '" "' ' dt 

where v*,.. = dv\\/dt = Vk,, / ^ = {k\\/k\\) is the group velocity of Bloch electrons with 
kinetic energy e^,, = hvpk\\ and the Fermi velocity vp ~ 10 8 cm/s . The Newton's second 
law requires that hd\i\\jdt = F^.. = — e (Ey + v^., x Bo) with — e, Ey and Bo being the 
electron charge, in-plane electrical and static magnetic fields, respectively. E||(r||, t) refers 
to the electric component of the total electromagnetic field, which should be determined by 
the self-consistent field equation (see next subsection below). The Markovian collision term 
introduced in Eq. ffTOl is 



d/(r||,k||,t) 
dt 



W£° [l-/(r,,,k„,t)]-VVg Hti /(r|[,k„,t) 



(out) 



(11) 



coll 
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where w£ n) and w£° ut) are the scattering-in and scattering-out rates for electrons in the 



two-dimensional kn-state. Based on the Newton's second law, we rewrite Eq. ([10]) into the 
form of 



df(r h k||, t) 
dt 



-v fc|l ■ V r „/(r||, k| h t) 



H 



Vk„/(r||, k,|, t) 



df(r h k h t) 



;i2) 



coll 



The zeroth-order moment of the Boltzmann equation in Eq. (1121) is found from 



2_ v 9/(r| h k||, *) 
A^ 



dt 



jYl v *n/( r ll> k b *) 



1 
h 



^Z) F *i- Vk »/( r ii' k ii'*) 



.4 



(13) 



k n k n 
which leads to the following electron number conservation equation, after the inter-valley 

scattering is ignored, 



dp(r h t) 
dt 



-Vr„ -j||(r||,t) 



(14) 



where A is the sample area, p(v\\, t) = (2 /A) J2f( r b k h ^ s ^ ne electron sheet number 

k " 

density (per area) and. j y (r y , t) = (2/A) ^2 v^. /(ry, ky, t) is the electron surface number cur- 

k 

rent density (per length). Equation (I14p allows us to determine the the spatial distribution 
of p(ry, t) at each time t. 

In order to simplify the first-order moment of the Boltzmann equation, we introduce the 
momentum-relaxation time approximation. Under this approximation, we write 



0/(rii,k[|,t) 



dt 



f( r b k ll> *) ~ /o( £ fcip T : A*o) 



(15) 



coll T l 

where /o(£fc||, ^\ A*o) = { ex P[( £ feii — Po)/^b^] + l} 1 is the Fermi-Dirac function for thermal- 
equilibrium electrons and T\ is the average momentum-relaxation time for electrons. In 
principle, T\ can be microscopically calculated based on Sjj>(t) introduced in Eq. (jl]) in the 
previous section for fixed applied bias and temperature as well as device parameters. In 
addition, we introduce the force-balance equation, which yields 
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(16) 



This leads to, to the leading order of a weak B field with vpB^jE^ <C 1, 



1 ~ eTl 

k|| « - — 



Em + u F 



Ey x B 



(17) 



Employing the result in Eq. (I15p . we arrive at the first-order moment of the Boltzmann 
equation 



-ti 



2 d 



v ^/( r ih k n^) 



4E v *ii (F fc|| • V k|| ) /(r,,, k„, 



A 



(18) 



Approximating /(ry, ky, t) on the right-hand-side of Eq. ffTSl by fo(ek v T, /iq), we get 



j||( r lh t ) + r i 



9j||(r||,*) 



~x}2 Wk \\ ( Vfc n ■ Vr n) /o( £ feii' T ' ^o) 



7-1 



2 ^ / \ df (e h , T, /i ) 



X] /o(£fc||, T, /i ) 



+ 



vperx 



Ey + Z/ F 



Ey x Bp 



2 >^ df (e kv T, no) 
A^ ~ 



9ei 



It is straight forward to show for each valley that 



dei 



(19) 



f (s k[[ , T, /i ) = A^bT) 2 Qi(r)) « p(ry, t) 



2 ^ a/o(efc||, T, /i ) 

-4 2^ h: = ~ N c ( k BT) Qo(v) 



de k 



(20) 
(21) 



where rj = /z //cbT, N c = \/{nh v F \ and the dimensionless function 
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x n dx 
3(35-77) _)_ 1 



(22) 



The local chemical potential /i, ( r || ? £) can De calculated from Eq. fl20|) if both T(r||, t) and 
p(ry, t) are given. Based on the results in Eqs. (120!) and (I2T1) . we finally get the generalized 
drift- diffusion equation 



0j||(r|h*) j||( r ll^) 



9i 



7i 



-Nn — 



En + z/ F 



En x B 



£11 



ObT) Qo(v) + - V r|| [(AfeT) 2 Q^rj)] j . (23) 

Equation ( 1231) enables us to determine the spatial distribution of the electron number 
current density j||(r||, t) at each time t. The Einstein relation can be obtained by set- 
ting <9p(r||, t)/dt = and <9j||(r||, t)/dt = in Eqs. ( !T4|) and (1231 for a steady state, i.e. 
V r|| ■ j||(r||, t) = 0, which relates the diffusion current to the external electric field E||(r||, t). 
In addition, by setting T(ry, t) = T L as a constant, Equations ({TBI and (1231) constitute 
the basic hydrodynamic equations for p(ry, £) and t). Although the particle number 

conservation is enforced in this way, the energy of the system is not conserved in general. 

The second-order moment of the Boltzmann equation is formally written as 



2^ <9/(r||, k| h t) 
k» 



dt 



-V P1 , • 



'2 

^E £k w v ^ii/( r ii' k u> *) 



1 
T, 



^E^ii F fc|, • V k|l /(rn, k,|, t) 



-■;(in) n _ t(^.. 1,., +W _ ^ c-. 1/\;(° ut ) 



A 



J^^W^f^^t). (24) 
k n " k n 

The results in Eq. (I24D can be simplified if we introduce the energy-relaxation time r 2 through 



e[T(r h t)]-e(T L ) 

T 2 



(25) 



where is the lattice temperature, and r 2 can be evaluated using the calculated non- 
equilibrium part of electron distribution g'j(t) as well as <Sjj>(t) in Eq. (j3J) for fixed applied 
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bias field, temperature and device parameters. This leads Eq. (1241) to the following electron 
power loss equation 



de[T(r h t)} e[T(r\\,t)}-e(T L ) 
- = V r|| ■ S||(r|j, t) H 



Of 

+ e 



En + vp 



Ej x B c 



72 

'j||( r M) 



(26) 



where the second and last terms on the right- hand- side of this equation corresponds to ther- 
mal energy exchange with lattice and Joule heating, e[T(ry, t)} = (2/ A) J2 £ ku f( r b k lh * s 

k n 

the average kinetic energy of electrons per area, and S||(r||, t) = (2/04.) J2 £ k^k ll /( r ||? k y, t) 

k u 

is the electron surface energy flux per length. It is easy to show that 



e[T(r h t)] « ^ J>*>ll /o(e*„ , T > *) = N c (k B T) 3 Q 2 (rj) 



(27) 

ky 

By substituting Eq. (|2"7|) into Eq. (l2"B"j) . this lets us find the spatial distribution of the electron 
temperature T(ry, t) at each time t. 

To simplify the third-order moment of the Boltzmann equation, we still employ the 
momentum- relaxation time approximation in Eq. ( II 5p . This leads to 



2 , d 

£fcn v fcn / ( r lh k lh *) + r i 



A 



dt 



-riV r 



2 x - 

~4 £fe n Vfe n Vfc ii ^( r H' k ll' *) 



A 



h 



~J £k W Vfe ll Ffe ll ' Vk ll^( r ll' k H' *) 



(28) 



In addition, approximating /(ry, ky, i) on the right-hand-side of Eq. (|28|) by /o (£&,,, T, /io), 
we get 



aSy(ry,0 S||(r||,t) 



at 



-N, 



En + vp 



Ey x Bp 



(W 2 Qi(r/) + ^V r|| [(A; B T) 3 Q 2 (v)]\ 



(29) 



From Eq. (|29|) . we are able to calculate the spatial distribution of the electron energy flux 
Sy(ry, t) at each time t, from which the electron thermal conductivity can be obtained. 
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Now let us summarize our findings in this subsection by presenting the complete set of 
equations describing the moments of the Boltzmann equation. Since Dirac electrons posses 
non-parabolic energy dispersion, the relevant moments are written as the powers of both 
energy and group velocity, i.e. 

M nl>n2 (r||, t) = 4EK,) nl K) n2 /(r||, k,|, t) , 
k n 

where M 0j o has the meaning of electron density, M 10 is the sheet current density, M 0) i is 
the average electron kinetic energy, and Mi i i is the electron surface energy flux. The above 
moments are augmented with the electron temperature T(r», t). Based on two calculated 
moments Mo,o and M^o, the momentum (or the group velocity) can be determined from the 
force-balance equation, i.e. ky = fa/ftjF^AMo^o; M 10 ), where the momentum relaxation- 
time fa) approximation has been employed. In general, the total force F^ , including the 
resistive force from electron scattering, must be found from the Maxwell equations. This 
point is further elucidated in the next paragraph. Note that the simultaneous solution of 
the Boltzmann-Maxwell equations is known as the Vlasov- Maxwell equation—, which can 
be used to elevate the relaxation-time approximation including the long-range Coulomb 
interaction. Our approach in this paper still retains the essence of the Maxwell equations 
but simplifies their introduction into the Boltzmann equation by means of the force-balance 
equation. As a result of this, it leaves us with the momentum and energy fa) relaxation 
times in the formalism. This partial simplification of the Vlasov- Maxwell equations gives 
rise to the correct treatment of electron transport driven by a terahertz optical field. 

The system of equations for the moments is given by 



M 0> o = -V r|| .Mi, , (30) 
M li0 = -^-iV c ^|F k|| (M 0i o;M li o)(A;BT) Q (r)) + - & V r|| [(k B T) 2 Q 1 (t 7 )]| , (31) 

M ,i = -V r|| • Mi i — _ eF ( M q) . M (32) 

r 2 

Mi,i = - 4e |iV c F k|| (M ,o; M 1)0 ) (k B T) 2 Q 1 ( v ) + lv r|| M 0il | , (33) 

where Mo,i = N^ksT) 2, Qzirj). As mentioned before, these equations are intertwined with 
the Maxwell equations if one aims to look for the self-consistent response of Dirac plasma 
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to an incident optical field. Another approach would be studying the ideal magneto- 
hydrodynamics of the plasma, which can be applied to the case when magnetic lines are 
frozen into an electron plasma. Formally, this corresponds to a vanishing total force 
(Ft,, = 0) and to a decoupling of the moments equations. Once those moments are obtained, 
it would provide the electromagnetic field inside the Dirac plasma via Maxwell equations as 
described in the next subsection. 

B. Self-Consistent Field Equation 

The Maxwell equations for the transverse magnetic component B(r, t) = B||(r, t) + 
Bj_(r, t) with B(rii, t) = B(rii,z = 0, t), as well as for the electric component E(r, t) = 
(r, t) + E_l(i-, t) with E(ry, t) = E(r||, z = 0, t), are given by 

V r ■ B(r, t) = , V r • [e r (r) E(r, t)\ = , (34) 



<9B(r, t) _ , <9E(r, t) 



at ' v ' x E(r - 4) • m 



e r (r 



2 

' V r x B(r, t) , (35) 



where e r (r) is the relative dielectric constant of embedded host materials including the gate 
oxide material and induced optical polarization field. The calculations of these four equations 
can be performed by using the Delaunay-Voronoi surface integration scheme.— The total 
electromagnetic fields, E(r, t) and B(r, t), are coupled to the moments of the Boltzmann 
equation through the following boundary conditions for {E z , By} at the two-dimensional 
graphene sheet [z = 0) 

e r (r||, + ) E z {r h + , t) - e r (r||, 0") E z (y h 0", t) 

[n s (r h t) + N ion -p(r h t)] , (36) 



e 



B y (r||, + , t) - B v (r h 0", t) = -efi oJx (r h t) , (37) 

B x (r||, + , t) - B x {r h 0", t) = e^j y {r h t) , (38) 
where A^ on is the ion sheet density, and the total charge neurility requires that 
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1 

A 



1 



d 2 r llP (v h t)=Ni 



ion 



(39) 



Here, the Maxwell equations must be solved self-consistently with the Boltzmann moment 
equations in the previous subsection.— In addition, we have continuity conditions at the 
boundary z = 



In Eq. ( l36i) . n s (r», t), which produces a space-charge field, is the induced surface charge 
density by a gate voltage Vc{t). For the graphene transistor structure, we also require that 

f dx E x (t\\, t) = Vvs(t) and f dzE z (r, t) = Vq(£), where we assume that the conduction 
o o 

channel is in the x direction with a channel length Lq and Vbs(^) represents the applied 
source-drain ac voltage. Also, L G represents the gate electrode depth. 

IV. CONCLUDING REMARKS 

In conclusion, we have found there is a minimum electron mobility for a graphene nanorib- 
bon just before a threshold for an applied electric field when entering into the nonlinear trans- 
port regime, which is attributed to the gradual build-up of a frictional force from phonon 
scattering by the applied field. We also predict a field-induced mobility enhancement right 
after this threshold value, which is regarded as a consequence of initially-heated electrons in 
high energy states with a larger group velocity in an elastic-scattering dominated graphene 
nanoribbon. Moreover, we have discovered that this mobility enhancement reaches a maxi- 
mum in the nonlinear transport regime as a combined result of an upper limit for the carrier 
group velocity in a nanoribbon and a field-induced dramatically-increased phonon scattering 
in the system. Finally, the threshold field can be pushed upward and the magnitude in the 
mobility enhancement can be reduced simultaneously by a small correlation length for the 
line-edge roughness in the high-field limit due to the occupation of high-energy states by 
field-induced electron heating. 



E x , y (r h + , t)=E x , y {r h 0", t) , 



(40) 



H z {r h + , t) = H z (r h 0", t) . 



(41) 



1(3 



Additionally, we have formulated a self-consistent device simulation for graphene transis- 
tors in the presence of field modulation within the terahertz frequency range. This involves 
moment equations from the Boltzmann equation up to the third order for fast carrier dynam- 
ics, as well as full wave electromagnetics coupled to the Boltzmann equation for describing 
both temporal and spatial dependence of the total field including the induced optical polar- 
ization field. 

When the electron density is increased in a graphene nanoribbon, multi-subband trans- 
port occurs, the field-induced mobility enhancement is expected to be reduced, and the effect 
of electron-electron scattering needs to be included. When the lattice temperature becomes 
high, on the other hand, neither the optical phonon nor inter-valley scattering should be 
neglected. As the width of a nanoribbon is modulated, a periodic potential along the ribbon 
will form, leading to a graphene nanoribbon super-lattice with additional mini-band gap 
opening at the Brillouin zone boundaries. The results of the current research is expected 
to be very useful for our understanding and design of high-power and high-speed graphene 
nanoribbon emitters and detectors in the terahertz frequency range. 
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FIG. 1: (Color online) (a) Calculated electron mobilities fi e as a function of applied electric field 
J-q at T = 10 K (solid squares on blue curve) and T = 6K (solid circles on red curve); (b) Electron 
drift velocities v& as a function of temperature T at J-q = 2kV/cm (solid squares on blue curve) 
and J-q = lkV/cm (solid circles on red curve). 
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FIG. 2: (Color online) (a) /i c as a function of J-q at T = 10 K with 70 = 1.0 x 10 13 s _1 (solid squares 
on blue curve) and 70 = 1.0 x 10 14 s^ 1 (solid circles on red curve); (b) Vd as a function of T with 
J-q = 2kV/cm for 70 = 1.0 x 10 13 s _1 (solid squares on blue curve) and 70 = 1.0 x 10 14 s~ 1 (solid 
circles on red curve). 



20 



0.45 

w 

> 0.40 

CM 

E 

o 

CO 

o 

H. 0.35 



0.30 



W = 40 A 
T = 10 K 


i i i i i 

y = 1 x 10 1 V 





^ • 




200 A 


- ^ A o = 


50A ^ „ 

y ^ 


(a) 

m m ■ — ■ — ■ — ■ H 


' n„ = 1 x 10 5 cm 1 

1D 



CO 



0.0 0.5 1.0 1.5 

F n ( kV / cm ) 



2.0 



1.0 
0.9 
0.8 
0.7 
0.6 
0.5 
0.4 



A = 200 A 
A = 50 A 



(b) 



W = 40 A 



n = 1 x 10 5 cm" 1 

1D 

y = 1 X10 13 S" 1 

F = 2 kV / cm 



10 



T (K) 



FIG. 3: (Color online) (a) fi c as a function of J-q at T = 10 K with Ao = 200 A (solid squares on 
blue curve) and Ao = 50 A (solid circles on red curve); (b) as a function of T with J-q = 2kV/cm 
for Aq = 200 A (solid squares on blue curve) and Ao = 50 A (solid circles on red curve) . 
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